There is room for goals to evolve based on feedback and scientific needs.
Load package
Pre-work was to install R, Rstudio, and pcvr with dependencies.
library(pcvr) # or devtools::load_all() if you are editinglibrary(brms) # for bayesian modelslibrary(data.table) # for freadlibrary(ggplot2) # for plottinglibrary(patchwork) # to arrange ggplots
Why Longitudinal Modeling?
plantCV allows for user friendly high throughput image based phenotyping
Resulting data follows individuals over time, which changes our statistical needs.
Why Longitudinal Modeling? 2
Longitudinal Data is:
Autocorrelated
Often non-linear
Heteroskedastic
Why Longitudinal Modeling? 3
Why Longitudinal Modeling? 4
Why Longitudinal Modeling? 5
Why Bayesian Modeling?
Bayesian modeling allows us to account for all these problems via a more flexible interface than frequentist methods.
Bayesian modeling also allows for non-linear, probability driven hypothesis testing.
Why Bayesian Modeling? 2
In a Bayesian context we flip “random” and “fixed” elements.
Fixed
Random
Interpretation
Frequentist
True Effect
Data
If the True Effect is 0 then there is an \(\alpha\cdot100\)% chance of estimating an effect of this size or more.
Bayesian
Data
True Effect
Given the estimated effect from our data there is a P probability of the True Effect being a difference of at least X
Why Bayesian Modeling? 3
Before moving on we should note that each group in your data will have parameters fit to it. If you have many groups then fitting models to only a few groups at a time is likely to make your life easier. Larger models will fit, but they can take a very long time.
It’s generally been easier for the authors to fit 50 models each with 4 groups in them then compare across models than to fit 1 model with 200 groups in it.
Supported Growth Models
There are 6 main growth models supported in pcvr.
3 are asymptotic, 3 are non-asympototic.
Supported Growth Models
Supported Growth Models
There are an additional 3 sigmoidal models based on the Extreme Value Distribution. Those are Weibull, Frechet, and Gumbel. The authors generally prefer Gompertz to these options but for your data it is possible that these could be a better fit.
Supported Growth Models 2
There are also two double sigmoid curves intended for use with recovery experiments.
Supported Growth Models 3
Generalized Additive Models (GAMs) are supported but discouraged due to poor interpretability.
Supported Growth Models 4
Intercept only models are supported, although they are only intended for use in segmented models or to represent homoskedasticity as a sub-model.
Supported Growth Models 5
Segmented growth models are also supported using “model1 + model2” syntax in both growthSim and growthSS.
For now we will focus on single models but there will be examples of segmented models later on.
Supported Growth Models 6
Finally, the “decay” keyword can be used to specify a decay model instead. For now we will focus on growth models as these constitute the large majority of models used in plant phenotyping.
Survival Models
Survival models can also be specified using the “survival” keyword. Using the brms backend these models use the weibull distribution by default but “binomial” family can also be specified by model = "survival binomial".
For details please see the growthSS documentation.
growthSS
Any of the aforementioned models can be specified easily using growthSS.
growthSS - form
With a model specified a rough formula is required to parse your data to fit the model.
id group time y
1 id_1 a 1 0.004058413
2 id_1 a 2 0.022060411
3 id_1 a 3 0.091828381
4 id_1 a 4 0.305290892
5 id_1 a 5 0.839871042
6 id_1 a 6 1.969862435
growthSS - form 3
Generally it makes sense to visually check that your formula covers your experimental design.
Note that it is fine for id to be duplicated between groups, but not within groups
growthSS - sigma
Recall the heteroskedasticity problem, shown again with our simulated data:
growthSS - sigma
There are lots of ways to model a trend like that we see for sigma.
pcvr allows any of the model options in growthSS to also be applied to model variance.
growthSS - Intercept sigma
Pros
Cons
Faster model fitting
Very inaccurate intervals at early timepoints
growthSS - Linear sigma
Pros
Cons
Models still fit quickly
Variance tends to increase non-linearly
Easy testing on variance model
growthSS - Linear sigma
growthSS - Gompertz sigma
Pros
Cons
Models fit much faster than splines
Slightly slower than linear sub-models
Variance is often asymptotic
Requires priors on sigma model
Easy testing on variance model
Note that these traits are broadly true of logistic and monomolecular sub models as well.
growthSS - Gompertz sigma
growthSS - Spline sigma
Pros
Cons
Very flexible and accurate model for sigma
Significantly slower than other options
Fewer priors
Splines can be a black-box
growthSS - Spline sigma
growthSS - other sigma models
You can always add a new sigma formula if something else fits your needs better.
growthSS - Distributional Models
Here we have limited the examples to talk about sigma, a parameter of the Student T distribution that our model belongs to. Written another way we might say all the previous methods are modeling:
Y ~ T(mu ~ main growth formula, sigma ~ sigma formula, nu ~ 1)
growthSS - Distributional Models 2
In pcvr the Student T family is the default for these models, but other distributions are supported through "distribution: model" syntax.
In general this is only for special cases where the Gaussian/T does not capture some important quality of the data. An obvious example could be leaf counts, which might be modelled as "poisson: monomolecular", for instance.
growthSS - Distributional Models 3
You can specify sigma as a list of formulas to model different distributional parameters separately. Most of the time this is overkill and adds unnecessary complexity, but the option exists for certain cases such as ZINB where modeling mean, shape, and zero inflation per group may make sense.
For details on supported families and their parameterization try running growthSS(..., sigma=NULL,...) and examining the priors, or checking ?brmfamily.
growthSS - priors
Bayesian statistics combine prior distributions and collected data to form a posterior distribution.
Luckily, in the growth model context it is pretty easy to set “good priors”.
growthSS - priors
“Good priors” are generally mildly informative, but not very strong.
They provide some well vetted evidence, but do not overpower the data.
growthSS - priors
For our setting we know growth is positive and we should have basic impressions of what sizes are possible.
At the “weakest” side of these priors we at least know growth is positive and the camera only can measure some finite space.
growthSS - priors 2
Default priors in growthSS are log-normal
\(\text{log}~N(\mu, 0.25)\)
This has the benefit of giving a long right tail and strictly positive values while only requiring us to provide \(\mu\).
growthSS - priors 3
We can see what those log-normal distributions look like with plotPrior.
Our final call to growthSS will look like this for our sample data.
ss<-growthSS(model="gompertz", form = y~time|id/group,sigma="gompertz", df=simdf,start =list("A"=130, "B"=10, "C"=0.5,"sigmaA"=20, "sigmaB"=10, "sigmaC"=0.25))
fitGrowth
Now that we have the components for our model from growthSS we can fit the model with fitGrowth.
fitGrowth 2
This will call Stan outside of R to run Markov Chain Monte Carlo (MCMC) to get draws from the posterior distributions. We can control how Stan runs with additional arguments to fitGrowth, although the only required argument is the output from growthSS.
fitGrowth 2
Here we specify our ss argument to be the output from growthSS and tell the model to use 4 cores so that the chains run entirely in parallel, but the rest of this model is using defaults.
fit <-fitGrowth(ss = ss, cores =4, iter =2000, chains =4, backend ="cmdstanr")
fitGrowth 2
Note that there are lots of arguments that can be passed to brms::brm via fitGrowth.
One that can be very helpful for fitting complex models is the control argument, where we can control the sampler’s behavior.
fitGrowth 2
adapt_delta and tree_depth are both used to reduce the number of “divergent transitions” which are times that the sampler has some departure from the True path and which can compromise the results.
fit <-fitGrowth(ss, cores =4, iter =2000, chains =4, backend ="cmdstanr",control =list(adapt_delta =0.999, max_treedepth =20) )
fitGrowth 3
fitGrowth returns a brmsfit object, see ?brmsfit and methods(class="brmsfit") for general information.
Within pcvr there are several functions for visualizing these objects.
growthPlot
growthPlot can be used to plot credible intervals of your model.
growthPlot(fit, form=ss$pcvrForm, df=ss$df)
Model Visualization 2
These plots can show one of the benefits of an asymptotic sub model well.
In pcvr some of what Kruschke recommends can be accessed from a model/list of models using the barg function, see ?barg for details on what that entails.
Using brms directly
These functions are all to help use common growth models more easily.
The choices in pcvr are a small subset of what is possible with brms, which itself is more limited than Stan.
It can be more work to try new options in brms or Stan, but if you have a situation not well represented by the existing models then it may be necessary.
Resources
If you run into a novel situation please reach out and we will try to come up with a solution and add it to pcvr if possible.
Good ways to reach out are the help-datascience slack channel and pcvr github repository.
